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Abstract 

The rational covariance extension problem to determine a rational spectral density given a finite 
number of covariance lags can be seen as a matrix completion problem to construct an infinite- 
dimensional positive-definite Toeplitz matrix the north-west corner of which is given. The circulant 
rational covariance extension problem considered in this paper is a modification of this problem to 
partial stochastic realization of reciprocal and periodic stationary process, which are better represented 
on the discrete unit circle 1j2n rather than on the discrete real line Z. The corresponding matrix 
completion problem then amounts to completing a finite-dimensional Toeplitz matrix that is circulant. 
Another important motivation for this problem is that it provides a natural approximation, involving 
only computations based on the fast Fourier transform, for the ordinary rational covariance extension 
problem, potentially leading to an efficient numerical procedure for the latter. The circulant rational 
covariance extension problem is an inverse problem with infinitely many solutions in general, each 
corresponding to a bilateral ARMA representation of the underlying periodic (reciprocal) process. In 
this paper we present a complete smooth parameterization of all solutions and convex optimization 
procedures for determining them. A procedure to determine which solution that best matches additional 
data in the form of logarithmic moments is also presented. 
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I. Introduction 

THE RATIONAL covariance extension problem or the partial stochastic realization problem 
has been studied in various degrees of detail in a long series of papers Pi-Ell, ifTTl . [|2TT|. 
Il22l . [|30l . [|39l . In a formulation suitable for this paper it can be stated as follows. Given a 
sequence (co, c%, . . . , c n ) of numbers, with c real and the rest possibly complex, such that the 
Toeplitz matrix 

C Ci c 2 ■ ■ • c n 
C\ C Ci • • • C n _i 
C~2 Ci C • • • C n _ 2 



(1) 



_C n C n _i C n _2 • • • Co 

is positive definite, find an infinite extension c n+ i, c n+3 , c n+3 , . . . such that, with c_fc 
k — 1, 2, . . . , the series expansion 



;6\ 



co 
k=— 00 



Cfc, 



(2) 



converges to a positive spectral density for all 9 G [—it, it) which takes the rational form 

P(z) 



Q(z) 



where P and Q are symmetric trigonometric polynomial of the form 



E 



Pke 



-ik9 



p-k = Pk, 



(3) 



(4) 



of degree n in the case of Q or at most n in the case of P. In [1211 . [|22l it was shown that 
there exists a Q for each assignment of P and in [OQ it was finally proved that this assignment 
is unique and smooth, yielding a complete parameterization suitable for tuning. Consequently, 
the rational covariance extension problem reduces to a trigonometric moment problem, where, 
for each P, the remaining problem is to determine a unique Q such that 

,,P{e %e ) d6 



Jk8- 



Cfc, k 



0,1,2, 



n. 



(5) 



Q{e id ) 2tt 

In 01, BH a convex optimization procedure to determine these Q was introduced, a result that 
has then been generalized in several directions 0-0, ll9l. ifTOl. IfTTl. E51, ll24t 
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The rational covariance extension problem can be seen as a matrix completion problem to 
construct an infinite-dimensional positive-definite Toeplitz matrix with T n in its north-west 
corner, which moreover satisfies the rationality constraint ©. There is an large literature on band 
extension of positive-definite Toeplitz matrices, starting with the work of Dym and Gohberg |[T6l 
and surveyed in the books ll25l . iflDl . dealing with the maximum-entropy solution corresponding 
to P = 1 . 

The circulant rational covariance extension problem, considered in this paper, is a modification 
of this problem to partial stochastic realization of periodic stationary process, which, as we shall 
explain in detail below, are better represented on the discrete unit circle Z 2 n (the integers 
modulo 2N) than on the the discrete real line Z. The corresponding matrix completion problem 
then amounts to completing a finite-dimensional Toeplitz matrix that is circulant. An important 
motivation for this problem is that its solution is a natural approximation of the solution to the 
ordinary rational covariance extension problem that turns out to involve only computations based 
on the fast Fourier transform and seems to lead to an efficient numerical procedure. 

Circulant matrices are Toeplitz matrices with a special circulant structure 







m 


m v 


m v -i ■ 


■ nil 






mi 


m 


m v 


■ m 2 


Circ{m , m 1 , m 2 , . . 


.,m„} = 


m 2 


nil 


m • 


■ m 3 






m u 


m u _i 


m v - 2 ■ 


■ m _ 



(6) 



where the columns (or, equivalently, rows) are shifted cyclically, and where m ,nii, . . . ,m u 
here are taken to be complex numbers. In the circulant rational covariance extension problem 
we consider Hermitian circulant matrices 



M := Circ{m , nii,m 2 , . . . , m N , m N _i, . . . , m 2 , nil}. 



(7) 



Hermitian circulant matrices appear naturally in the context of periodic stationary stochastic 
processes. To see this consider a zero-mean stationary process {y(t)}, in general complex- valued, 
defined on a finite interval [-N+1, N] of the integer line Z and extended to all of Z as a periodic 
stationary process with period 2N; i.e., such that y(t + 2kN) = y(t) almost surely. Processes of 
this kind can naturally be defined on the group Z 2N of the integers with arithmetics modulo 2N, 
and in this setting stationarity can be seen as propagation in time of random variables under the 
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action of a (finite) unitary group. We shall write the string {y(— N + 1), . . . ,y(0), . . . ,y(N)} 
as a 2iV-dimensional column vector y and only consider stationary process of full rank, whose 
covariance matrix 

£:=E{yy*} (8) 

is positive definite, where * denotes transpose conjugate. Let K{y(t) \ y(s), s ^ t} be the wide 
sense conditional mean of y(t) given all {y(s), s ^ t}. Then the error process 

d(t):=y(t)-E{y(t)\y(s),s^t} (9) 

is orthogonal to all random variables {y(s), s ^ t}; i.e., E{y(t) d(s)} = a 2 5 ts , t,s e Z 2A r, 
where 5 is the Kronecker function and a 2 is a positive number, or, equivalently, 

E{yd*} = a 2 I, (10) 

where I denotes the 2N x 2N identity matrix. Interpreting © in the mod 2N arithmetics of 
%2N, y admits a linear representation of the form F y = d, where F is a 2N x 2N circulant 
matrix with ones on the main diagonal. Following Masani [37 ; ], d is called the (unnormalized) 
conjugate process of y. Therefore, setting e := -^d and A := -^F, we see that a full-rank 
stationary periodic process admits a normalized representation 

Ay = e, E{ey*} = I, (11) 

where A is Hermitian and circulant. Since AE {yy*} = E {ey*} = I, A is also positive definite 
and the covariance matrix © is given by 

S = A" 1 , (12) 
which is circulant, since the inverse of a circulant matrix is itself circulant. Therefore, if 

c k := E {y(t + k)yjt)}, k = 0, 1, 2, . . . , N, (13) 
S is precisely the Hermitian circulant matrix 

S = Circ{c , ci, c 2 , . . . , c N , c N ^, . . . , c 2 , c x }. (14) 

In fact, a stationary process y is full-rank periodic in Z 2 at, if and only if S is a Hermitian 
positive definite circulant matrix |fTT|. 
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We are now in a position to state the main problem of this paper. Supposing that only the 
covariance lags c , Ci, . . . , c n are available for n < N, how do we complete the matrix (fl4l) with 
the entries c n+ i, c n+2 , . . . , Cjy so that it is circulant and the covariance matrix © of a stationary 
periodic process y with a spectral density of the rational form ©. We would like to parametrize 
the set of all solutions to this problem. 

This can be seen as a generalization of modeling of reciprocal processes about which there 
is a large and important literature flU, ESI, (2*0, 03, 03-051 . A first step in this direction 
was taken in ifTTTl . where the circulant matrix A in (fT2l is required to be banded of order n; i.e., 



A = Circ{a , a x , . . . , a n , 0, . . . , 0, a n , a n _i, . . . , a x }. 
For example, a banded matrix of order n = 2 takes the form 



(15) 



A 



a 


a i 


«2 










a 2 


ai 


a i 


«o 


ai 


a 2 










«2 




a i 


a 


fli 


«2 















ai 


a 


«1 


a 2 















«2 


«1 


«o 


Oi 


a-2 












«2 


a i 


a 


a i 


a L 


«2 










«2 


ai 


a 



(16) 



In this case y admits a bilateral AR representation 

n 

^2 a kV(t - k) = e(t) , 

k=—n 

for all t G 1<2N, which can also be written 

A(()y(t) = e(t) 

in terms of the symbol 



ci-k — °>k 



M0 = E a ^~ k ' 



at — au 



(17) 



(18) 



(19) 



k=—n 



of A, where ( represents the forward shift on Z 2 n- These models are representations of stationary 
reciprocal processes of order n; see |fTT|. where they are determined by solving a maximum- 
entropy problem. 
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In this paper we show that for each choice of banded circulant matrix P of order at most n, 
there is an unique banded circulant matrix Q of order n such that 

S = Q X P. (20) 

If the corresponding symbols are P(C) and Q((), respectively, then, by (PTTI) and (PT21) . such a 
solution corresponds to a bilateral ARMA representation 

Q{()y{t) = P(C)e(t), (21) 

or, equivalently, 

n n 

g fc y(t - k) = ^2 Pk e ( f - k )- ( 22 ) 

k=—n k=~n 

We have therefore a complete parameterization of such representations, and hence of the com- 
pletions of S, in terms of the P. 

In Section [TTl we review basic facts about circulant matrices and harmonic analysis on Z 2 n 
and set up notations. The main results on the complete parameterization of the circulant rational 
covariance extension problem are presented in Section Hn] where we also consider the circulant 
rational covariance extension problem as an approximation procedure for the ordinary rational 
covariance extension problem. In Section [TV] following flU, JH, ifTTll . lfT8ll . [|24ll . we show how 
the parameter P can be determined from logarithmic moments computed from data. 

II. Preliminaries 

Most of the harmonic analysis of stationary processes on Z carries over naturally, provided 
the Fourier transform is understood as a mapping from functions defined on Z 2 tv onto complex- 
valued functions on the unit circle of the complex plane, regularly sampled at intervals of length 
A := n/N. We shall call this object the discrete unit circle and denote it by T 2 at. This Fourier 
map is usually called the discrete Fourier transform (DFT). Next we shall review some pertinent 
facts and at the same time set up notations. 

A. Harmonic analysis on Z 2A r 

Let Ci : = e* A be the primitive 2iV-th root of unity; i.e., A = tt/N, and define the discrete 
variable £ taking the 2N values ( k = = e lAk ; k = —N + 1, . . . ,0, . . . , N running counter- 
clockwise on the discrete unit circle T 2 n- In particular, we have = (complex conjugate). 
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The discrete Fourier transform J maps a finite signal g = {g k ; k = —N + 1, . . . , N}, into a 
sequence of complex numbers 



-N + 1, -N + 2, . . . , N. 



(23) 



N 

G(G):= E ' • / 

It is well-known that the signal g can be recovered from its DFT G by the formula 

- A 

9k= E CjGiQ^r, k = -N + l,-N + 2,...,N, (24) 

j=-iV+l 

where = ^ plays the role of a uniform discrete measure with total mass one on the discrete 
unit circle T 2A r. In the sequel it will be useful to write (|24|) as an integral 

Qk = e ike G{e ike )dv{6), k = -N + 1, -N + 2, . . . , N, 



(25) 



where v is a step function with steps ^ at each ( k ; i.e. 



N 



-jV+1 



2iV 



(26) 



It easy to check that, if F is the DFT of {ft}, 

N N 

E /*& = ^ E F ((kMU) = / F(e^)G(e^)*<M0). 

fe=-JV+l k=-N+l ^ n 

This is Plancherel's Theorem for DFT. 

It is sometimes convenient to write the discrete Fourier transform (T2~3l in the matrix form 



(27) 



g = Fg, 

where g := (G(C- N+ i), G(C- N+2 ), • • • , G(( N ))\ g 
nonsingular 2N x 2N Vandermonde matrix 



(28) 

(g-N+i, 9-N+2, • • • , 9n) T and F is the 



C 



iV-1 r N-2 
-N+l S-/Y+1 



>/V-2 

So So 



5JV S/Y 



a- TV 



-iV 



-iV 



(29) 



Likewise, it follows from (1241) that 



F*£ 

2N 8 ' 



(30) 



i.e., ? 1 corresponds to n^F*. Consequently, FF* = 2NI, and hence F 1 



1 ^F*and(F*)- 1 



J-F 
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B. Circulant matrices 

A Hermitian circulant matrix 

M := Circ{m , mi, m 2 , . . . , tun, Tn>N-i, ■ ■ ■ , ^2, mi}- 
can be represented in the form 

TV 

k=-N+l 

where S is the nonsingular 2N x 2N cyclic shift matrix 

1 ... 

1 ... 

1 ... 

S := 

1 

1 j 

which is itself a circulant matrix with symbol S(Q = (. Clearly S 2N = S° = I, and 



(3D 



(32) 



(33) 



ik+2N 



*2N-k 



s- k = (S 



k\T 



Consequently, 



SMS* = M, 



(34) 



(35) 



and the condition (1351 ) is both necessary and sufficient for M to be circulant. (Clearly, what is 
said in this section holds for circulant matrices in general, but in this paper we are only interested 
in the Hermitian ones.) 

As before setting g := (g- N+ i, g~ N +2, • • • , 9n) t , we have 



[Sg]fc = gk+i, k e Z 2N - 
In view of (T23T) . it then follows that C3 r (g)(^) = 3 r (Sg)(C), from which we have 

j(M g )(0 = M(CMg)(C), 

where the pseudo-polynomial 

N 



(36) 



(37) 



(38) 



fc=-JV+l 
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is called the symbol of the circulant matrix M. 

An important property of circulant matrices is that they are diagonalized by the discrete Fourier 
transform. More precisely, it follows from (|37l) that 

M = ^F*diag(M(C-Ar +1 ), . . . , M(C-i), M(Co), M(G), . . . , M(( N ))F, (39) 

and consequently the inverse M _1 is 

M- 1 = ^F*diag(M(C-^v+i)- 1 , . . . , M(Ca,)- 1 )F. (40) 

Since 

S = ^F*diag(C- J v + i,...,G)F and S* = ^F*diag(C^ +1 , . . . , Q l )F, 

we have 

SM _1 S* = M _1 . 

Consequently, M _1 is also a circulant matrix with symbol M(£) -1 . In general, in view of the 
circulant property (|32l) and (|34l) . quotients of symbols are themselves pseudo-polynomials of 
degree at most N and hence symbols. The coefficients of the corresponding pseudo-polynomial 
M(£) _1 can be determined by Lagrange interpolation. More generally, if A and B are circulant 
matrices of the same dimension with symbols A(Q and B(() respectively, then AB and A + B 
are circulant matrices with symbols A(QB(() and A(()+B(Q, respectively. In fact, the circulant 
matrices of a fixed dimension form an algebra and the DFT is an algebra homomorphism of 
the set of circulant matrices onto the pseudo-polynomials of degree at most iV in the variable 

C e T 27V . 

C. Spectral representation of periodic stationary stochastic processes 

Let {y(t)} be a zero-mean stationary process defined on a finite interval [— N + 1, N] of 
the integer line Z and extended to all of Z as a periodic stationary process with period 2N. 
Let C-jv+i, C-n+2, ■ ■ ■ , cat be the covariance lags Ck := E {y(t + k)y(t)}, and define its discrete 
Fourier transformation 

N 

E c ^i^ j = ~N + l,...,N, (41) 

k=-N+l 
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which is a positive real-valued function of £. Then, as seen from (1241) and (|25l) , 

- a r 

c *= E ^(0)3-= / e^$(e ifce )^), fe = -iV + l,...,iV. (42) 
The function $ is the spectral densitiy of the process 7/. In fact, let 

N 

V((k):= E y(t)C fc _t , fc = -iV+l,...,iV, (43) 

t=-JV+l 

be the discrete Fourier transformation of the process y. The random variables (l43~l) turn out to 
be uncorrelated, and 

^{v((k)v(CtY} = HCk)Skt. (44) 

This can be seen by a straight-forward calculation noting that 

1 - 

2^ E = (45) 

t=-JV+l 

Then, we obtain a spectral representation of y analogous to the usual one, namely 

V(t)= E C^(a)^=/ e^dyCC), (46) 

where 

#(0) := y{e i6 )dv{e). (47) 

It is interesting to note that 

= 2^^{y(a)y(a)*}, fc = -N + 1, . . . , N, (48) 
obtained from (|44b has a similar form as the periodogram 

HCk) = ^y((k)y((k)*, k = -N + i,...,N, (49) 

widely used in statistics. 
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III. The complete solution to the circulant rational covariance extension 

PROBLEM 

Given n < N and cq, c%, . . . , c n with a positive definite Toeplitz matrix (OQ), find a spectral 
density $ satisfying the moment conditions 

e lk %(e w )dv = c k , k = 0, 1, 2, . . . , n. (50) 

Note that this moment condition can also be written as an underdetermined system of linear 
equations 

1 N 

^ E = fc = 0,l,2,...,n, (51) 

j=-N+l 

in the variables = J = — ^ + !>•••? — 1> 0, 1, • • • , N, where the coefficient matrix is a 

Vandermonde matrix of full rank. Note that it is consistent with (|50l) to define negative moments 
by setting c_fc = Ck, so that the pseudo-polynomial 

n 

C(C) = E c k C k (52) 

fc=— n 

is the symbol of a banded Hermitian circulant matrix 

C = Circ{c , ci, . . . , c„,, 0, . . . , 0, c n , c„_i, . . . , c x } (53) 

of order n. We would like to find a rational extension c n+ i, c n+ 2, . . . , cat replacing the zeros in 
C to obtain a Hermitian circulant matrix 

S := Circ{c , c x , c 2 , . . . , c N , cjv-i, • • • , c 2 , (54) 

that is positive definite. In terms of stationary periodic processes this is the covarance matrix 
([8]). We proceed to solve this in terms of the symbols, and then interpret the results in terms of 
matrices. 

A. Circulant rational covariance extension in terms of symbols 

Let ^3 be the finite-dimensional space of symmetric trigonometric polynomials ©, and define 
to be the positive cone 

<p + = {p e *p | P(e id ) > for all 9 e [-7T, tt]}. 
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Moreover, let £ + be the dual cone of all c = (c , c x , . . . , c n ) such that 

n 

(C, P) := c kPk > for all Pg^, (55) 

k=—n 

where the notation (C, P) is motivated by the fact that, by dTTT) . 

(C,P)= / C(e ie )P(e i0 )*dis. (56) 



It can be shown that c G <£+ if and only if T„ > 0. In fact, if a(z) = a z n + ■ ■ ■ + a„_iz + a n 
is a polynomial spectral factor of P{z), i.e., a(z)a(z)* = P(z), then it is easy to see that 

(C,P) = a*T n a, (57) 

where T n is the Toeplitz matrix (Q} of c , ci, . . . , c n and a = (a , ai, . . . , a n ) T - 
Next, define the cone 

qj + (JV) = {Pe«p|p(c*)>o fc = -iv + i,-iv + 2,...,iv}. (58) 

Clearly, ( P+(A r ) D $P+(2iV) D $P+(4iV) D ■ ■ • D and the corresponding dual cones satisfy 

<t+(N) c £+(2iV) c £+(4iV) C • • • C €+. (59) 

Theorem 1: Let c G C + (N). Then, for each P 6 ^3 + (iV), there is a unique Q G ^3 + (iV) such 
that 

* = I (60) 

satisfies the moment conditions (|50l ). 

For the proof, which is given in the appendix, we need to consider a dual pair of optimization 
problems. First consider the primal problem to maximize the generalized entropy 

fir 

I P ($) = / P(e i0 )\og<$>(e ie )dv (61) 



subject to the moment conditions (1501) . The corresponding Lagrangian is then given by 

n / pit \ 
L($,Q) = I P ($) + J2 ^ [c k - / e m ^ e )dA 

k=-n ^ , '- 7r ' 

= r P(e* e ) log $(e ie )^ + (C, Q) - r g(e i9 )$(e^)rfz/, (62) 

J —TV J— 7T 
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where q , q±, . . . , q n are Lagrange multipliers, and where Q is defined as in © with q_ k = q k . 



Since the dual functional sup$ Q) is finite only if Q G *}3 + (iV), we may restrict the Lagrange 



multipliers to that set. Therefore, for each Q G ty + (N), consider the directional derivative 
which equals zero for all variations <5$ if and only if 

Q 

Inserting this into (l62l) we obtain 

P(e id ) [\ogP(e ie )-l] dv, 

where 

J P (Q) = (C, Q) - f P(e w ) \ogQ(e w )du (63) 

and the last term is constant. Hence we may take (l63l) as the dual functional. 

It will be shown below that Jp is strictly convex, so a stationary point in if it exists, 
would have to be a unique minimizer of Jp. For k = 1, 2, . . . , n, we write q k = x k + as a 
sum of real and imaginary parts and define the partial differential operators 

«t^— ] (64a) 

(64b) 



dq k 2 \dx k dy k 
d 1 ( d d 



dq k 2 \dx k dy k 
in the standard way; see, e.g., Il27l p. 1]. It is immediately seen that 



dqk _ dq k 

— = and — = 0. (65) 



From this, we readily obtain 



Setting (|66l ) equal to zero yields the moment conditions (|50l ). Then the proof of the following 
theorem follows directly from Theorem [Q 

Theorem 2: Let Let c e <£+(A r ) and P G ^ + (iV). Then the problem to maximize (|6"TT) subject 
to the moment conditions (|50l ) has a unique solution, namely (|60l) , where Q is the unique optimal 
solution of the problem to minimize (|63l) over all Q G ^3 + (iV). 
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From (|66l we have the Hessian 



r-= M = 0,l,...,n, (67) 



J_„~ Q(e ie ) 2 

which is Hermitian and positive definite, showing that Jp is strictly convex. 

Next, we establish that the solution to the circulant rational covariance extension problem 
depends smoothly on the parameters c and P. To this end, for each fixed P E y$ + (N), we 
define the moment map F p componentwise given by 

r p(p i6 ) 

K(Q) = / e^-^-J-du, k = 0, 1, . . . , n (68) 



Q{e ie 

and, for each c E €+, the map G c sending P E <P+(iV) to Q E Q+(N) := (F p )-\€ + (N)). 
The proof of the following theorem is given in the Appendix. 
Theorem 3: The maps F p and G c are homeomorphisms. 

In particular, we have established a complete smooth parameterization of all solutions Q to 
the circulant rational covariance extension problem in terms of P E ^3 + (iV). 

Finally, as a prerequisite for Theorem [7] in the next section, we show that the map F p can 
be continuously extended to the boundary < ^ + (A r ), as can be seen from the following extension, 
proved in the Appendix, of the family of dual solutions. 



Theorem 4: Let c E <t + (N), and let ( p + (A r ) denote the closure of *}3 + (iV). Then, for each 



P E *P+(iV) \ {0}, the dual problem to minimize d63]) over Q E *P+(iV) \ {0} has a unique 
minimizer Q, and P/Q satisfies the moment conditions (|50l ). 

B. Circulant rational covariance extension in terms of matrices 

Next we reformulate the optimization problems in terms of circulant matrices. To this end, 
we define the circulant matrix 

S = ^F*diag($(C^+i), . . . , HCn))F (69) 
with symbol (|60l ) and the banded numerator matrix 

P = ^F*diag(P(C-iv + i), • • • , P({n))F (70) 

of degree at most n with symbol (HJ). Since $(Cfc) > for all k and log$(C) is analytic in the 
neighborhood of each $(00 > 0, by the spectral mapping theorem [15, p. 557] the eigenvalues 
of log S are just the real numbers log ^(Cfc), k = —N + 1, . . . , N, and hence 

logS = ^F*diag(log$(C- JV+ i), • • • ,log$(6v))F. (71) 
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Consequently, the primal functional (|6TT) may be written 



PIT 1 N 

/ P{e ie )\og<b(e w )dv = — J2 ^(0)log$(C 

J A— AT I 1 



7V+1 

^tr{Plo g 5]} (72) 



and the moment conditions (1501) as 



^tr{S fc S} = c fc , fc = 0,l,...,n, (73) 



or, equivalently, as 

E„ T SE n = T n , where E n = ' 



(74) 





Consequently the primal problem amounts to maximizing tr{P log £} over all Hermitian, positive 
definite 2N x 2N matrices subject to (1731) or (|74|) . For the special case P = 1 this reduces to the 
primal problem presented in [[TT|. except that in 0TJ there is an extra condition insuring that £ 
is circulant. However, it was shown in [12J that this condition is automatically satisfied and is 
therefore not needed. 

In the same way, by (l56l) the dual functional (|63l) can be written 

C(e id )Q(e i9 )dv - / P(e w ) log Q(e i9 )du 

J-n (75) 

= ^tr{CQ}-^tr{PlogQ}, 

where 

Q = ^F*diag(g(C-iv + i),...,g(Civ))F (76) 

and C is the banded circulant matrix (1531) formed from Co, ci, . . . , c n . 

Consequently, given c G £+(iV), it follows from Theorem Q] that, for each Hermitian, positive 
circulant matrix P that is banded of degree at most n, there is a unique £ given by 

£ = Q X P, (77) 

where Q is the unique solution of the problem to minimize 

Jp(q) = ^ tr { c Q} - ^ tr i p lQ g Q} ( 7g ) 

over all q := (q , q x , . . . , q n ) such that the Hermitian, circulant matrix 

Q = Circ{g , ?i, • • • , q~n, 0, . . . , 0, q n , qn-!, . . . , a x } 
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is positive definite. For the maximum-entropy solution corresponding to P = I this reduces to 
an optimization problem that is different from the one presented in ifTTI . 
Since 

n 

Q = qkS q - k = 

k=— n 

we have 

= S and = S fc , 

for fc = 0, 1, . . . , to, and therefore 

if = 2^< slc > " M tt ^ p<i ^ = <* - 4 tt { s ' PQ ">' (79) 

where we have used the fact that 



tr 



{s fc c> = c M sk ~ j } = 2Nc k 



j=-n 



as tr{S fe } 7^ only for S° = I. Setting (|79b equal to zero yields the moment conditions. Likewise, 

<92jp - 1 tr{S fe PCT 2 S- £ } = -Ltr{S k - e PQ- 2 }, (80) 



2N L 2AT 
showing that the Hessian is a Toeplitz matrix. This is the matrix version of (|67l ). 

In IfTTI it was observed that the condition that the Toeplitz matrix T n , defined by ©, is 
positive definite is a necessary, but not a sufficient, condition for feasibility of the circulant 
banded covariance extension problem. This can now be understood in the more general setting 
of moment problems discussed above. In fact, the Toeplitz condition T„ > is equivalent to 
c £ £+. whereas, by Theorem [2l c £ €+(N) is required for feasibility. Since C + (N) C <£+, it 
follows that the Toeplitz condition cannot be sufficient in general. However, as proved in [[TT||. 
feasibility is achieved for a sufficiently large N. This can also be seen from the following result. 

Proposition 5: The feasibility set £+(N) — > £ + as N — > oo. In particular, for any c £ <£+, 
there is an A" such that c £ £+(N) for N > iVo- 

Proof: As N — > oo, the set {Q; j = —N + 1, . . . , N} becomes dense on the unit circle, 
and therefore ^ + (A^) — > *$ + . Consequently, € + (N) — > <£ + , and the convergence is monotone in 
the sense of (l59l) . Therefore, since £ + is an open set, there is an N such that any c £ £ + will 
sooner or later end up in £+(AT) and remain there as A" increases. ■ 
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C. Some computational considerations 

For each choice of P, the Hessian of Jp can be computed explicitly in terms of Q as the 
Toeplitz matrix 



(81) 





h 


h 


h 2 ■ 


K 




h 


ho 


hi ■ 


h n -i 


H(q) = 


h 2 


k 


h ■ 


' h n _ 2 




K 


hn-l 


h n ~2 ■ 


h 



where q = (g , qi, ■ ■ • , q n ) are the coefficients in the pseudo-polynomial Q and 



N 



— T c k - 

j=-N+l ^ S J' 



(82) 



as can be seen from (|67l) or (f8Qb . Therefore Newton's method can be used to find the unique 
minimizer of the dual problem. The gradient (|66|) at the point q is (c — c(q)), where 



c(q) = — V C k ^^- 

j=-N+l ^^J^ 

and consequently a Newton step amounts to solving the Toeplitz system 

H(q W ) (q (fc+1) -q W )=c-c(qW). 
Clearly, H(q) and c can be computed by the discrete Fourier transform. 



(83) 



(84) 



D. An approximation procedure for the ordinary rational covariance extension problem 

Given a c E £ + and a P E the ordinary rational covariance extension problem amounts 
to finding the unique Q E satisfying the moment conditions 

c * = lj aP <mt k=0A >-- n <85) 

We would like to approximate the solution Q of this problem by the unique solution Qn of the 
circulant rational covariance extension problem 

pf e iB\ 



Ck 



ike 



where dujsr is the measure (T2~6l) corresponding to N 



di^N, k — 0, 1, . . . , n, 



(86) 
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A N) I r^I^l^-, I, 0.1 „ (87) 



Theorem 6: Let c G £+ and a P £ ^} + . Moreover, for any iV > iVo, where N is defined 
as in Proposition [5j let Q N be the unique solution of (f86l) . and let C} be the unique solution of 
(1851 Then Q as JV -»■ oo. 

Proof: Let F : — > C + be the map sending Q to c as in in (IB3T ); i.e., c = F(Q). Given 
Qjv> define c* 7 ^ := F(Q N ) with components 

S N ) _ / ^ifce 

~Q N {e i0 )2ir'- 

for each N > N . Since (l86l) is a Riemann sum converging to (f87l) as in the measure dfjy 
tends to oo but Qm is kept fixed, there is for each e > an N\ > N such that ||c( w ) — c|| < e for 
all N > Ni. Consequently, since F(Q N ) = c^- N \ F(Q) = c and the map F is a diffeomorphism 
[|71 Theorem 1.3], Qn — > Q as in iV — > oo. ■ 

IV. Determining P from logarithmic moments 

We have shown that the solutions of the circulant rational covariance extension problem are 
completely parameterized in a smooth manner by the numerator trigonometric polynomials P e 
*P + (iV), or, equivalently, by their corresponding banded circulant matrices P. Next, we show 
how P can be determined from the logarithmic moments 

e ^log$(e ie )dz/, fc = l,2,...,n. (88) 

In the setting of the classical trigonometric moment problem such moments are known as cepstral 
coefficients, and in speech processing, for example, they are estimated for design purposes. 

Now consider the problem of finding the spectral density $, or, equivalently, the circulant 
matrix S, that maximizes the entropy gain 

I($) = y log$(e*Vz/=^trlog£ (89) 

subject to the two sets of moment conditions (|50l) and (f88l . Such a problem was apparently first 
considered in the usual trigonometric moment setting in an unpublished technical report ||38l 
and then, independently and in a more elaborate form, in fl5j, J6j, ifTTll . 
Defining 

n 

M(C) = m ^" fc , (90) 

k=—n 
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where m_ k = m k , k = 1, 2, . . . , n, and m = 0, the Langrangian for this optimization problem 
can be written 

n / pit \ 

P, Q) = I($) + Y^q k [c k - e m ^(e ie )du 

k=-n ^ ' 
n / j-K \ 

-^2pk(m k - / e <fcfl ]og$(e*)di/J 

fe=-n V J ~* ' (91) 

= (C,Q) - (M,P) - ^ Q(e id )<S>(e w )dv 
+ / P(e id )\og$(e w )dv, 

J —TT 

where pi, . . . ,p n , g , q x , . . . , q n are Lagrange multipliers and p := 1, and P and Q are the 
corresponding trigonometric polynomials ©. For the dual functional (P, Q) h-> sup $ P, Q) 
to be finite, P and Q must obviously be restricted to the closure of the cone <p + (iV). Therefore, 
for each such choice of (P,Q), we have the directional derivative 

5L{$, P, Q- 5$) = J \^-QJ8$dv, (92) 



and hence a stationary point must satisfy 



• = (93) 



which inserted into (1911 ) yields 



where 



supL($,P,Q) = J(P,Q)-1, 



J(P, Q) = (C, Q) - (M, P) + P(e te ) log SJy*', (94) 

and where we have used the fact that J" Pdf = po = 1- Accordingly, we define the bounded 
subset 

(AT) := {P e ^ + (iV) I = I}- (95) 

of the cone *p + (Af). Note that J is convex, but not necessarily strictly convex unless P and Q 
are coprime, and that 

w* = Ck ~L e gbr k = 1 — n (96a) 

° ] r e ike log^J-du-m k , fc = l,...,n. (96b) 
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Consequently, if there exists a stationary point (P, Q) £ (N) x <p + (iV), (|93l will satisfy both 
the moment conditions (|5Q|) and (l88l) . 

A proof of the following theorem, which is a circulant version of Theorem 5.3 in will be 
given in the Appendix. 

Theorem 7: Suppose that c £ £+(-/V) and mi, . . . , m n are complex numbers. Then there exists 
a solution (P, Q) that minimizes J(P, Q) over all (P, Q) £ *#%{N) x ^p+(iV)» and, for any such 
solution 



satisfies the covariance moment conditions (|5Q|) . If, in addition, P £ *}3 + (iV), (|971) also satisfies 
the logarithmic moment conditions (f88l) and is an optimal solution of the primal problem to 
maximize the entropy gain (|89l) given (|50l) and ((881) . Then Q £ Vp + (N), and the solution is 
unique. In fact, J is strictly convex on *}3° (iV) x *}3 + (iV). 

Consequently, solving these optimization problems will always lead to a spectral density 
with the prescribed covariance lags c , ci, . . . , c„, provided c £ <£ + (iV). However, we have not 
prescribed any condition on the logarithmic moments mi, ... , m n , as such a condition is hard to 
find and would depend on c. If the moments c , Ci, . . . , c n and mi, ... , m n come from the same 
theoretical spectral density, the optimal solution (|97l) will also match the cepstral coefficients. In 
practice, however, Co, ci, . . . , c n and mi, ... , m n will be estimated from different data sets, so 
there is no guarantee that P does not end up on the boundary of *}3 + (iV) without satisfying the 
logarithmic moment conditions. Then the problem needs to be regularized, leading to adjusted 
values of mi, . . . , m n consistent with the covariances c , c\, . . . , c n . 

We shall use a regularization term proposed by P. Enqvist ifTTl in the context of the usual 
rational covariance extension problem. More precisely, we consider the regularized dual problem 
to find a pair (P, Q) £ ^° + (N) x *P+(JV) minimizing 



P 



(97) 



Q 




(98) 



for some A > 0, or in circulant matrix form 



Ja(PQ) = ^tr{CQ} " ^tr{MP} + -i-trjPlogPQ- 1 } - -^trllogP}. (99) 
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This functional will take an infinite value for P G <9^3 + (iV), since then P(Ck) — for some k, 
and hence the minimum will be in the interior. Then 

X -' ' k6 log-±J-du-m k -e k = 0, k = l,...,n, (100) 



dpk & Q(e ie ) 

at the minimum, where 



and hence the moments (pOl) and (|88l) are matched provided one adjusts the logarithmic moments 
mi,m 2 , . . . ,m„ to mi + £i,m 2 + e 2 , • • • , m « + £n> the latter of which are consistent with 
Co, Ci, . . . , c n . Modifying the analysis in ifTTl p. 188 - 196] to the present setting it is easy 
to see that (l99l is a monotonically nonincreasing function of A, and that the solution tends as 
A — > oo to a (P, Q) where P = 1, i.e., the maximum entropy solution. 
Computing the Hessian of Ja, we notice that 



dq k dq t 2iV. = ^ +i ^ Q(C 



is the same as (1671) . Moreover, 

N 



dq k dp e 2N^ +l * 3 Q((j 



Af N 



d3x = J_ V f fc ~' 1 +— V f fc ~' A (102c) 

Since J is strictly convex (Theorem[7]), then so is Ja, so the Hessian is positive definite. Newton's 
method can then be used as in Section IIII-CI to determine the unique minimizer. 

V. Conclusions 

In this paper we have presented a complete parameterization of all solutions to the circulant 
covariance extension problem. We have shown that determining these solutions involves only 
computations based on the fast Fourier transform, potentially leading to efficient numerical 
procedures. This also provides a natural approximation for the ordinary rational covariance 
extension problem. 

The circulant rational covariance extension problem is an inverse problem with infinitely many 
solutions in general, but by matching additional data in the form of logarithmic moments a unique 
solution can be determined. 
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For many applications it will be important to generalize these results to the multivariable case. 
This should be straight-forward, but we have chosen to consider only the scalar case in this paper 
in order to keep notations reasonably simple and not blur the picture. 

Appendix 

A. Proof of Theorem [7] 

Consider the moment map F p : ^ + {N) — > € + (N) defined by (1681) for an arbitrary P E 
y$ + (N). This is a continuous map between connected spaces of the same (finite) dimension. 
Therefore, if we can prove that F p is injective and proper - i.e., for any compact K C £ + (N) 
the inverse image (F p ) -1 (fQ is compact - then, by Theorem 2.6 in [8], it is a homeomorphism, 
implying in particular that the system of moment equations F P (Q) = c has a unique solution 
in ?p + (JV). 

Lemma 8: The moment map F p : ty + (N) — » £+(iV) is injective. 

Proof: From (1631) we have the gradient (l66l) and the Hessian (|67l) . which is positive definite. 
Therefore, Jp is strictly convex, and any stationary point is a solution to the moment equations 
(l50l) . which must be a unique if it exists. Hence F p is injective. ■ 
It remains to show that there exists a solution to the moment equations (l50l) . 
Lemma 9: Suppose the Toeplitz matrix T n is positive definite; i.e., c £ <£ + . Then, for any 
compact K C € + (N), the inverse image (F P )^ 1 (K) is bounded. 

Proof: Suppose Q satisfies the moment equations F P (Q) = c for some c E £ + (TV). Then 

n /"7T 

(C,Q) = c ^ = / P{z W )dv =: «, 

where k is a constant. Now, let a(z) = a z n + ■ ■ ■ + a n -\Z + a n be the stable polynomial spectral 
factor of Q(z), i.e., a(z)a(z)* = Q(z). Then 

K={C,Q) = a*T nS L, 

where T„ is the Toeplitz matrix of c and a ) '. If c is restricted to the compact 

subset K E £ + , the eigenvalues of T n are bounded away from zero. Hence T n > el for some 
e > 0, and consequently 

a 2 < -a*T n a= -. 

e e 

Consequently, ||q||, where q = (g , ?i, • • • , q n ), is a l so bounded, and hence so is (F p )~ l (K). ■ 
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Lemma 10: The moment map F p : ty + (N) — > € + (N) is proper. 

Proof: Let K be a compact subset of £ + (N), and let c( fe ) be a sequence in K converging 
to c G if. Since (F p ) _1 (if) is bounded (Lemma [9]), there is a convergent sequence in 
the preimage of the sequence c( fc ) converging to some limit Q. We want to show that Q e 
(F P ) -1 (.K'). The only way this can fail is that Q belongs to d^ + {N), the boundary of yS+(N). 
We observe that 



"7T p2 

(C {k \P) = I ^Tjrdv, 



, Q ik) 



and consequently, since P E *$ + (N), 



j=-N+l VlSjJ 

which requires that Q(Q) ^ for all j. However, Q can only belong to dty + (N) if some Q(Q) 
equals zero. Hence Q ^ dty + (N), as required. ■ 
This concludes the proof of Theorem [IJ 

B. Proof of Theorem \3\ 

We have already proven above that F p is a homeomorphism. It remains to prove that G c is. 
For this we need two more lemmas. 

Lemma 11: For each fixed c G £+(N), the map G c : 0.+ (N) — > *}3 + (iV) is injective. 

Proof: Suppose that G C (Q 1 ) = G C (Q 2 ) = P for some Qi,Q 2 E Q+(N). We want to show 
that Qi = Q 2 . To this end, since 

j —-a 

we have 



-dv = 0, fc = 0,l, 

(Qi-Q 2 ) 2 i' J .._ ^ A, (i) „ ( 2)l r.«(«i-ft) p 



where 



fc=— ?i 

and, consequently, Qi(Cy) = ^MCj) f° r ai l J» as claimed. 

Lemma 12: For each fixed c e £ + (iV), the map G c : 0.+ (N) — > ^} + (7V) is proper. 
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Proof: The proof follows the same pattern as that of Lemma [101 Hence, let K be a compact 
subset of *}3 + (iV), and let P^ be a sequence in K converging to P G K. Since (G c ) _1 (i^) C 
H + (iV) is bounded (Lemma [9]), there is a convergent sequence in the preimage of the 
sequence P^ converging to some limit Q. In order to ensure that Q G (G )" 1 ^), we must 
demonstrate that Q G' d*$ + (N). To this end, note that 



and consequently, since c G £+(N) and P £ K C %$ + (N), 

jr S^ = (c,p)>o, 

j=-iV+l $(0) 

Since P((j) > for all j, this requires that Q(Q) > for all j. Hence Q £ d^ + (N), as 
required. ■ 
The map G c is a continuous map between connected spaces of the same dimension n + 1. 
Noting that (1671) is positive definite, the continuity follows from the inverse function theorem 
applied to the equation F P (Q) = c. Then, since G c is injective and proper, it follows from 
Theorem 2.6 in [8], that it is a homeomorphism. 

C. Proof of Theorem |4] 

We follow the lines of the proof of the Main Lemma in ||7l p. 10]. Let (Pi) be a sequence 



in ( p + (A r ) converging to P G y$ + (N) \ {0}. Then there is a positive constant K such that 
Pt((j) < K for i = 1, 2, 3, . . . and j = -N + 1,...,N. For each £, let Q e be the unique 
minimizer of 

J Pi (Q) = (C,Q) - r P l (e j£ )\ogQ(e j£ )dv 

J — TT 

as prescribed by Theorem [2l Then 

e* k9 ^^dv = c k , fc = l,2,...,n. (103) 

Now suppose that the sequence Qe is unbounded. Then there is a subsequence, which we shall 
also denote Qt, for which ||Q^||qo > 1 an d ||Q^||oo — > oo. For each such Q, there is an e > 
such that 

> e||Q||oo - i^log HQlloc. (104) 
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To see this, first note that, since T n > 0, (C, Q/||Q||oo) has a minimum e > on the compact 



set {Q e ¥+(N) | HQIU = 1}, we have (C,Q) > e||Q||oo. Then 

JpM > e||g|U - f p,\og (j^j dv - Kl °z IIQIloo, 

where the second term is nonnegative and can be deleted. 

Next, let Q E *#+{N) be arbitrary. Then, by optimality, l Pl (Q) > $p t (Qe)- Since $p t (Q) -> 
Jp(Q) as t — > oo, there is a positive constant L such that 

L> Jp £ (Q) > I Pe (Qi), £=1,2,3,..., 

which together with (11041) yields 

£>£||Q£||oo-^log||Q^|U ^=1,2,3,.... (105) 

Then, comparing linear and logarithmic growth, we see that the sequence (Qt) is bounded, 
contrary to hypothesis. Consequently, there is a convergent subsequence (for convenience also 
indexed by £) such that Q t — > Q, and, since Pi — > P ^ 0, (11031) implies that Q ^ 0. Hence, 
setting $^ := P?/Q^ and $ := P/Q, — > and hence, taking limits in (11031) . we obtain 

e ike $(e ie )dv = c fc , fc = 1,2,..., ra, 
showing that Q is the required minimizer satisfying the moment conditions. 

D. Proof of Theorem [7| 

We begin by showing that the sublevel set JJ -1 (oo, r] is compact for each rGl. The sublevel 



set consists of those (P, Q) e ¥l( N ) x ¥+( N ) for which 

r > h{P,Q)+h{P), 

where 

h(P,Q) = (C,Q) - f P(e ie )\ogdetQ(e ie )du 

J — 7T 

/7T 
P(e i6 )logdetP(e ie )dv 
-IT 

Since *}3° (N) is a bounded set that is bounded away from zero, there is a positive constant K 
such that || P ||oo < K and a p 6 R such that J 2 (P) > p for all P e (iV). Hence, in view of 
the estimates leading to (11041) . 

r-p> Ji(P, Q) > ellQlU - K\og HQIU, 
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and therefore, comparing linear and logarithmic growth, it follows that the sublevel set JJ _1 (oo, r] 
is bounded. Since it is also closed, it is compact, as claimed. 

Since J thus has compact sublevel sets, there is a minimizer (P,Q). Then clearly Q is a 
minimizer of Jp, and hence, by Theorem H <f> := P/Q satisfies the moment conditions (|50l) . 
If P G *$?(N), then the minimizer must satisfy the stationarity condition dj/dpk — 0, k — 
1, 2, . . . , N, and hence, by (|96b| ), <f> also satisfies the logarithmic moment conditions (|88l ). Since 

I(P, Q) = L{P, Q) > L{P, Q) for all (P, Q), 

and L(P, Q) = I(P, Q) for all (P, Q) satisfying the moment conditions (l50l) . (P, Q) solves the 
primal problem. By Theorem [Q Q G P+(iV). 

It remains to prove that the optimal solution is unique if (P, Q) G *P+(iV) x *p + (AT). To this 
end, we form the directional derivative 

5S(P, Q; 6P, 5Q) = (C- PQ~\ 5Q) + (MPQ- 1 ) - M, 5P) 

and the second directional derivative 

5 2 J(P, Q; 5P, 5Q) = ((SP- PQ- l 5Q), p-\6P - PQr x 6Q)) > 

with equality if and only if 

5P - PQ- l 5Q = 0. 

Then, however, 



/7T /*7T 
PQ- 1 5Qdu= / 5Pdv = 0, 
■7T ^/ — 7T 



since the pseudo-polynomial <5P has no constant term, as P(0) = 1. Therefore, choosing SQ = 1, 
it follows from Theorem Q] that 



7T 

-1 



Co = / PQ- 1 dv = 0, 

^ — 7T 

which is a contradiction. Consequently, 

5 2 J(P, Q; 5P, 5Q) > 

for all 8P, 5Q; i.e., the Hessian of J is positive definite, and hence J is strictly convex. Therefore 
uniqueness follows. 
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